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Abstract. This article concerns two methods for reducing large systems of chem- 
ical kinetics equations, namely, the method of intrinsic low-dimensional manifolds 
(ILDMs) due to Maas and Pope [U. Maas and S. B. Pope, Combustion and Flame 88 
(1992) 239-264] and an iterative method due to Fraser [S. J. Fraser, J. Chem. Phys. 
88 (1988) 4732-4738] and further developed by Roussel and Fraser [M. R. Roussel 
and S. J. Fraser, J. Chem. Phys. 93 (1990) 1072-1081]. Both methods exploit the 
separation of fast and slow reaction time scales to find low-dimensional manifolds 
in the space of species concentrations where the long-term dynamics are played out. 
The analysis is carried out in the context of systems of ordinary differential equa- 
tions with multiple time scales and geometric singular perturbation theory (GSPT). 
A small parameter e measures the separation of time scales. The underlying assump- 
tion is that the system of equations has an asymptotically stable slow manifold 
in the limit as e J, 0. Then it follows from GSPT that there exists a slow manifold 
Me for all sufficiently small positive e, which is asymptotically close to Mq. 

It is shown that the ILDM method yields a low-dimensional manifold whose 
asymptotic expansion agrees with the asymptotic expansion of up to and in- 
cluding terms of C(e). At 0{e^), an error appears that is proportional to the local 
curvature of Mq] it vanishes if and only if the curvature is zero everywhere. 

The iterative method generates, term by term, the asymptotic expansion of the 
slow manifold ^Ae■ Starting from Mq., the ith. application of the algorithm yields 
the correct expansion coefficient at 0(e*), while leaving the lower-order coefficients 
at 0{1) through 0(e*~^) invariant. Thus, after £ applications, the expansion is 
accurate up to and including the terms of 0{e^). 

The analytical results are illustrated in two examples: a planar system from 
enzyme kinetics (Michaelis-Menten-Henri) and a model planar system due to Davis 
and Skodje. 
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1 Introduction and Summary of Results 



Many chemical reaction mechanisms in combustion |j73|, atmospheric 
science |6^, enzyme kinetics |12], and biochemistry involve large numbers 
of species, multiple chains of chemical reactions, and widely disparate time 
scales. A typical model of hydrocarbon combustion, for example, may well 
involve several hundred species, which participate in hundreds of reactions 
that proceed on time scales ranging from nanoseconds to minutes. The size 
and complexity of these mechanisms has stimulated the search for methods 
that reduce the number of species and chemical reactions but retain a desired 
degree of accuracy. Typically, these reduction methods select a small number 
of species, which are marked as reaction progress variables, and determine the 
concentrations of the remaining species as functions of the latter, either by 
table look-ups or by direct computation. The critical step in these methods 
is, of course, the definition of the reaction progress variables, which may be 
actual concentrations of selected species or combinations thereof. 



Research into reduction methods has increased dramatically over the past 
decade, and several methods have been proposed in the hterature and imple- 
mented in computer codes. We mention the quasi-steady state approxima- 
tion ||62|, |7l|, the partial-equilibrium approximation ||66|, methods based on 
details of the chemistry [BB, E8|, an iterative method O, p3|, the method of 



intrinsic low-dimensional manifolds [S^] , the computational singular pertur- 
bation method |T8|, |19|, |32|, a principal-component analysis, lump- 
ing techniques repro-mo deling |j70|, an inertial- manifold approach a 
dynamic dimension- reduction method |^, a saddle-point method a 
predictor-corrector method an optimization method and a global- 
eigenvalue method |65]. 



In this article, we focus on two reduction methods, namely, the intrinsic 
low- dimensional manifold (ILDM) method due to Maas and Pope j^B], ^ 



and the iterative method due to Eraser and further developed by Roussel 
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and Fraser |5^. Both methods have been developed for and extensively 
applied to problems with slow manifolds that attract nearby initial conditions. 
The long-time behavior of such systems is governed by the dynamics on the 
slow manifold, whose dimension is generally much less than that of the total 
composition space, resulting in a considerable reduction of complexity. 

Given the importance of slow manifolds, a central question for any reduc- 
tion method is: How accurately does it approximate a slow manifold? The 
present investigation answers this question for the ILDM method of Maas and 
Pope and the iterative method of Fraser and Roussel. 

In the ILDM method, the Jacobian of the vector field is partitioned at 
each point of phase space into a fast and a slow component, and bases for the 
corresponding subspaces are generated by means of a Schur decomposition. 
The intrinsic low-dimensional manifold is defined as the locus of points where 
the vector field lies entirely in the slow subspace and is an approximation of 
the slow manifold. The efficacy of the ILDM method is evident, for example, 
by the reduction achieved in the prototypical example of a CO-H2-O2-N2 com- 



bustion model 123) HI- Disregarding only the production of NO, the model 
comprises evolution equations for the enthalpy, pressure, and concentration 
of each of thirteen species, making for a fifteen-dimensional phase space, and 
a total of sixty-seven chemical reactions. With the proper choice of the re- 
action progress variable (CO2), a reduction to a one-dimensional ILDM can 



be achieved that retains a certain accuracy after an initial transient |3^] 



Reduction to a two-dimensional ILDM gives a better approximation, albeit at 
the expense of keeping track of two reaction progress variables and the stor- 
age of a correspondingly larger look-up table. Refinements, applications, and 
evaluations of the ILDM method against direct numerical simulations can be 
found in Refs. [||, |3|, |T0|, |0|, |T], |n[ |l 0, |7|, |, |5|, ^ |I|, |, |7|, |5|, |] 



The iterative method was inspired by the phase-space geometry of an 
enzyme-kinetics model involving a fast and a slow species, where the slow 
manifold is a curve in the phase plane. The method is derived formally from 
the invariance equation — an equation that is satisfied on any trajectory of the 
dynamical system and, in particular, on the slow manifold and extends nat- 
urally to multidimensional systems with (possibly) higher-dimensional slow 
manifolds. The procedure is explicit if the force field is linear in the fast 
variable, and implicit otherwise; hence, it generally requires the use of a 
nonlinear equation solver. The method has been developed further and ap- 
plied, in particular, to several problems of enzyme kinetics and metabolism in 



Refs. H, P, m, 0, B 55, m 



A natural framework for the analysis of these and similar reduction meth- 
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ods is provided by geometric singular perturbation theory (GSPT) |T^ , ^ 



P^ , |25[| . The presence of a fast and a slow time scale leads naturally to the 
introduction of a small positive parameter e measuring the ratio of the char- 
acteristic times. If, in the limit as e J, (infinite separation of time scales), the 
system of kinetics equations has a slow manifold, phase space and this 

manifold is asymptotically stable, then GSPT identifies a (usually nonunique) 
slow manifold A^^ for e sufficiently small positive and gives a complete ge- 
ometric and analytical description of all solutions in the vicinity of the slow 
manifold, including how trajectories approach the manifold. By comparing 
the asymptotics of the slow manifold M.e found by GSPT with the asymp- 
totics of the low-dimensional manifolds generated by the ILDM method and 
the iterative method we can evaluate the accuracy of these reduction methods 
for small values of e (finite but large separation of time scales). The evaluation 
leads to the following conclusions. 

ILDM Method, (i) The asymptotic expansion of the ILDM agrees with the 
asymptotic expansion of the slow manifold M./, up to and including the 
0{e) term, for all fast-slow systems. In general, however, the 0{e'^) 
terms differ. 

(ii) The error at 0{e'^) is proportional to the local curvature of the slow 
manifold A^o- K vanishes if and only if the curvature of M.q is zero 
everywhere. (The "if" part was observed previously in Ref. [37].) 



Iterative Method, (i) The iterative method, if started from M.o, generates 
term by term the asymptotic expansion of the slow manifold M.^- In 
particular, I applications of the iterative method generate an approxi- 
mation to the slow manifold M./, that is asymptotically correct up to and 
including the 0{e^) term, albeit with extraneous terms at 0{e^^^). 

(ii) The £th iteration leaves the terms at 0{1) through 0{e^~^) invari- 
ant. (This observation is important because the lower-order terms have 
already been determined correctly in the preceding iterations.) 



Remark. In Ref. ||5^, it is shown that the ILDM coincides with the slow 
manifold M.q in the limit of infinite separation of the fast and slow time scales 
(e = 0). I 



Remark. The slow manifold M.q can often be found analytically; otherwise, it 
can be obtained by one application of the iterative method to the steady-state 
approximation. The latter is readily found numerically; see Refs. [pTS] , |53| . I 

The present article is organized as follows. In Section 0, we review the 
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general framework of fast-slow systems of ordinary differential equations and 
recall the asymptotic expansion of the slow manifold. In Section ^, we define 
the ILDM and indicate briefly how it is computed. We present the asymptotic 
expansion of the ILDM for planar fast-slow systems (one fast and one slow 
variable) in Section ^ and for general fast-slow systems {n fast and m slow 
variables) in Section ^. The results are summarized in Corollary |5.1| . In 
Section we describe the iterative method of Fraser and Roussel. We discuss 
its asymptotics in Section |^. The results are summarized in Corollary |7.1| . 
We illustrate the analytical results with two planar examples, namely the 
Michaelis-Menten-Henri mechanism of enzyme kinetics (Section ^ and an 
example due to Davis and Skodje (Section ^). In Section |10| we remark on 
several generalizations and discuss some remaining issues. 



2 Fast-Slow Systems of ODEs 

We consider reaction mechanisms in homogeneous media, where the concentra- 
tions of the chemical species depend on time only. The concentrations evolve 
on two distinct and widely separated time scales. The slowly evolving con- 
centrations are the entries of the vector y, the remaining concentrations the 
entries of the vector z; the former has m components, the latter n [m, n > 1). 
The separation of time scales is measured by e, an arbitrarily small positive 
parameter. The limit e | corresponds to infinite separation. The reac- 
tion mechanism is thus modeled by a system of ordinary differential equations 
(ODEs), 

y' = 6f{y,z,e), (2.1) 
z' = g{y,z,e). (2.2) 

The unknowns y and z are functions of t with values in R"* and R'^, respec- 
tively; ' denotes differentiation with respect to t; and / and g are smooth 
functions with values in R"^ and R", respectively. We assume that / and g, 
as well as all their derivatives, are C(l) as e | 0. 

Remark. The system of Eqs. ( |2.1| )-( p^ ) is, of course, an idealization of the 
complex systems that occur in chemical kinetics. The model is adopted here 
because it is suitable for mathematical analysis. We claim, however, that it 
also captures the essential elements of any reaction mechanism whose long- 
term dynamics evolve on slow manifolds and offers a paradigm for the analysis 
of reduction methods. The validity of our conclusions extends therefore well 
beyond the idealized system of Eqs. (|2.1| )- (|2.2| ). We comment on the implica- 
tions for more realistic systems in Section |10|. I 
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The independent variable t is called the fast time because it defines the 
time scale on which the fast variables evolve, and the system of Eqs. (|2.1| )- 
( |2.2|) is labeled the fast system. While the fast time scale is appropriate for the 
study of the transient dynamics, the long-time dynamics are more naturally 
studied in terms of the slow time r = et. On the scale of r, the system of 
Eqs. (|2rT| ) -(p^ ) assumes the form 

y = f{y,z,e), (2.3) 
£z = g{y,z,e). (2.4) 

Here, ' denotes differentiation with respect to r. We refer to the system of 
Eqs. (p73|)-(p^) as the slow system. 



The fast system ( |2.1|) - (|2.2|) and the slow system (|2.3|) - (|2.4|) are, of course. 



equivalent as long as e > 0, but they approach different limits as e I — that 
is, as the separation of the fast and slow time scales becomes infinite. The fast 
system reduces to 

y' = 0, (2.5) 

z' = giy,z,0), (2.6) 

which is essentially a single equation for the fast variable z with y as a param- 
eter. The slow system, on the other hand, reduces to 

y = f{y,z,o), (2.7) 

= g{y,z,0). (2.8) 

The first equation describes the motion of the slow variable y, and the second 
equation is an algebraic constraint that forces the motion to take place on the 
zero set of g. 

Our focus is on systems for which the zero set of g is represented by 
the graph of a function. That is, we assume that there exists a single-valued 
function h^, which is defined on a compact domain K = [0, F]™ in R*", such 
that 

g{y,ho{y),0) = 0, y E K. (2.9) 
The zero set of g thus defines a manifold, Aio, in phase space. 

Mo = {{y, z) e R'"+" : z = h^iy), y e K}, (2.10) 

to which the motion of the reduced slow system is confined. 

Our analysis requires a second assumption that holds for many, though 
not all, of the systems in which reductions have been sought, namely, that 
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each point {y, ho{y)) on JHq is an asymptotically stable fixed point of Eq. (p.6|). 
The assumption guarantees that the eigenvalues of the matrix Dzg{y, ho{y), 0) 
all have negative real parts. 

Remark. The two assumptions are justified in most enzyme kinetics and some 
combustion and atmospheric chemistry problems. In certain more complex 
reaction mechanisms, however, they may need justification. Toward this end, 
we observe that, in those cases where reduction methods are expected to be 
effective, Hq can be found locally by the Implicit Function theorem (since the 
second assumption guarantees that the matrix {Dzg){y, hQ{y),0) is invertible 
for each y E K), and GSPT can be applied to each local portion. In the 
absence of singularities, these local functions can be pieced together to form a 
smooth global function over the entire domain under consideration. | 



Under the above conditions, standard asymptotic theory (see, for exam- 
ple, Refs. |3^, |l^ I, ^ |2^) guarantees that, when e is positive but 
arbitrarily small, there exists a slow manifold M.,, that is invariant under the 
dynamics of the system of Eqs. (|2.1|) -( |2r2D , has the same dimension as A^o? 
and lies near /Aq. All nearby solutions relax exponentially fast to A^^, and 
their long-term evolution is determined by an associated solution on the slow 
manifold itself. The manifold M.^ is usually not unique; there typically is a 
family of slow manifolds, all exponentially close {0{e~'^/^) for some c > 0). 



Theorem 2.1 (Fenichel, asymptotically stable slow manifolds). For any suf- 
ficiently small e, there is a function that is defined on K such that the 
graph 

M, = {{y,z):z = K{y). V & K} (2.11) 

is locally invariant under the dynamics of Eqs. \2.]\ )- (^^ ). The function 
admits an asymptotic expansion, 

K{y) = hoiy) + 6h^'\y) + e''h^^\y) + ■■■ as e j 0, (2.12) 

where the coefficients h^^^ : K R" are found successively from the equation 

{Dzg)h^'^ = |:(D/.«)/(^-i-) - Y: ^(Dig) J] (/.(-), . . . , Z.^)) 
i=o i=2^- \i\=e 

-E^E^Pi(^e^7)) E {h^'^\...,h('^^)-^{d',g){2.13) 
k=i j=i 3- \i\=e.-k ^■ 

for i = 1,2,..., with h'^'^^ = Hq. Here, the functions f and g and their deriva- 
tives are evaluated at {y, z = ho{y), 0), and it is understood that a sum is empty 
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when the lower bound exceeds the upper bound. In particular, h^^^ and h^'^^ are 
given by 



{D,g)h^'^ = {Dho)f-ge, (2.14) 



- liDlg) {h^'\ /^(i)) - {D,g,)h^'^ - \g,,. (2.15) 

Furthermore, G C'^{K) for any finite r, and the dynamics of the system of 
Eqs. ^2.1\ )-(^^) on M.^ are given by the reduced equation 

y = f{y,hM,e). (2.16) 



Proof. The theorem is a direct restatement of Theorem 2 in Ref. ||2^ for 



the special case in which A^o is asymptotically stable. It also follows directly 
from p5| , Theorem] and is a special case of the Fenichel theory |T^. The 



asymptotics of the slow manifold Ais are given explicitly, for example, in 



Refs. p|, pi. I 



Remark. In many instances — for example, in the Michaelis-Menten-Henri 
reaction mechanism discussed in Section ^ and various combustion problems — 
the reduced slow system y = f{y,ho{y),0) has an asymptotically stable fixed 
point at (2/05^0(2/0)), say. In such cases, the reaction scheme has a global 
attracting equilibrium. Under the hypotheses made above, the system of 
Eqs. ( p.l| )-( p^ has a fixed point at (2/o,e; ^ed/o.e))? and the slow manifold 
Ais is its weak stable manifold. | 

Remark. While we have used it here only for the case of attracting manifolds. 



the Fenichel theorem and Theorem 2 in Ref. p2[ hold for the more general 
case of fast-slow systems of ODEs for which the manifold A4o is normally 
hyperbolic — that is, where there can be both fast stable (exponentially con- 
tracting) and fast unstable (exponentially expanding) dynamics in the direc- 
tions transverse to A4o- In the more general case, the matrix {Dzg){y, holy), 0) 
has s eigenvalues with a negative real part and u eigenvalues with a positive real 
part, the fast variable z decomposes into a u-dimensional and a s-dimensional 
component with u + s = n, and the dynamics of all solutions near Ai^ are 
governed by the Fenichel normal form The asymptotics of Ate remains 
unchanged. I 



Remark. The articles of Tikhonov [^] and Levinson |2^ present the origi- 
nal theory of persistence of asymptotically stable manifolds (see also Ref. [^ ). 
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The theory of persistence of normally hyperbolic manifolds can be found in 
the monographs of Fenichel [|l^, [1^ and Hirsch, Pugh, and Shub |^ ; see also 



Ref. [72]. Other relevant references are [45| and |iT] for singularly perturbed 
systems of ODEs with asymptotically stable slow manifolds and and 
for singularly perturbed systems of ODEs with general normally hyperbolic 
slow manifolds. An introductory exposition of GSPT is given in Ref. . I 



Remark. A numerical procedure for finding asymptotically stable slow mani- 
folds in fast-slow systems, which is stable and highly accurate for small values 
of e, has been given by Nipp |H6[. | 



3 The ILDM Method of Maas and Pope 



The ILDM method starts from the slow system, Eqs. (|2.3| )- (|2.4| ), takes the 
local vector field F and the associated Jacobian J, and reduces the latter at 
each point to a fast and a slow component. The vector field F and its Jacobian 
J are 

fJ I], J=i Ri ). (3.1) 

\e 9 J \e Dy9 e Dz9 J 

where Dyf is the m x m matrix of partial derivatives dfi/dyj, D^f the m x n 
matrix of partial derivatives dfi/dzj, DyQ the n x m matrix of partial deriva- 
tives dgi/dyj, and Dz9 the n x n matrix of partial derivatives dgi/dzj. 

By assumption, the real part of each eigenvalue of J is negative. The sum 
of the eigenvalues is equal to the trace of J, which is 0{e~^) as £ | 0, and 
their product is equal to the determinant of J, which is 0{e~"') as £ | 0. The 
eigenvalues of J fall therefore into two groups: one group of m eigenvalues 
with 0{1) negative real parts and another group of n eigenvalues with 0{6^^) 
negative real parts. The eigenvectors associated with the first group span the 
slow subspace, those associated with the second group the fast subspace. The 
Maas and Pope algorithm defines the ILDM as the locus of all points {y, z) 
where the vector field F lies entirely in the slow subspace. 



The algorithm uses a Schur decomposition Section 6.3] of J, 

J = QNQ', (3.2) 

with Q unitary {QQ' = Q'Q = Im+n, ' denoting the transpose) and N upper 
triangular, 

Q = {QsQf), N=(^^^ ^^f). (3.3) 



/ 
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The dimensions of Qs and Qf are (m + n) x m and (m + n) x n, respectively; 
A^s is an m X m upper triangular matrix, Nf annxn upper triangular matrix, 
and Nsf an m x n full matrix. The eigenvalues of J appear on the diagonal 
of N in descending order of their real parts, from least negative at the (1, 1) 
position to most negative at the {m + n,m + n) position. This particular 
ordering is accomplished in Ref. |^ by means of a modification of Stewart's 
implementation of the Schur algorithm |^ and in Ref. ||3^ by means of a 



standard Schur decomposition followed by a sequence of Givens rotations [17, 
Section 5.1]. 

The first m Schur vectors — that is, the columns of Qs — form an orthogonal 
basis for the slow subspace, while the remaining n Schur vectors — the columns 
of Qf — form an orthogonal basis for the orthogonal complement of the slow 
subspace. The vector field F is entirely in the slow subspace if it is orthogonal 
to the orthogonal complement of the slow subspace — that is, if 

Q'fF = 0. (3.4) 

This equation defines the ILDM, the latter being an approximation of the slow 
manifold Aie- We analyze its asymptotics (as e | 0) in the following sections. 

Remark. The matrix Q'r corresponds to Qj^, the number n to nj, and the 
sum m + n to n in Ref. . I 



In the numerical implementation of the ILDM method for general, closed, 
adiabatic, and isobaric reaction mechanisms, the system of equations is closed 
by supplementing the ILDM equation, Eq. ( p.4|) , by a set of parameter equa- 
tions. The parameter equations fix the enthalpy, pressure, and element compo- 
sition. In addition, the reaction progress variables are treated as parameters. 
Each fixed set of parameters yields one point of the ILDM, and the entire 
ILDM is obtained by sweeping over the admissible set of parameter values. 
As noted in Ref. [^, the parameters can generally be chosen so the ILDM is 
at least defined piecewise, and, most important, the choice of the parameter 
equations does not influence the construction of the manifold. 

In Eqs. ( p.l|) -( p72D , the enthalpy, pressure, and conserved quantities have 
been neglected. In this case, the parameter equations fix the values of the slow 
variables y, and the ILDM is obtained by sweeping over all points y & K. 
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4 Asymptotics of the ILDM — Planar Case 



We first restrict our attention to planar fast-slow systems, Eqs. (|2.3|) - (|2.4|) 
witfi m = n = 1, for which the computations are relatively straightforward 
and the asymptotic analysis more transparent. We address the general case in 
Section |. 

In the planar case, the vector field F and its Jacobian J are 

F=( i ], j=( ^ V (4.1) 

\^ 9 J \ £ 9y £ 9z J 

The eigenvalues of J are 

Kf = \ {^~'9z + /,) ± ^J\{e-^9z + fyf-e-^ ify9z-fz9y), (4.2) 

where the upper (lower) sign is associated with (A/). Thus, 

Xs = fy-^ + 0{6), Xf = e-'g, + 0{l) as £ i 0. (4.3) 
9z 

The derivatives of / and g, which are evaluated at {y,z,e), are all 0{1) as 
e I 0. The (nonnormalized) slow eigenvector is 



As - £ ^9z 



(4.4) 



and there is a corresponding fast eigenvector vj. The vector Vg spans the 
slow subspace, Vf the fast subspace. The vectors Vg and vj are not necessar- 
ily orthogonal. To determine the points [y, z) in the phase plane where the 
vector field F lies entirely in the slow subspace, we work with the orthogonal 
complement of the slow subspace, which is spanned by the row vector 

= {e-^gy, e-^g, - A^). (4.5) 

The locus of all points in the phase plane where the vector field F is in the 
slow subspace coincides with the set of all points (y, z) where F is orthogonal 
to f ^ — that is, where 

fgy + g{e''gz-\s)=Q. (4.6) 
This equation defines the ILDM. 

Remark. A Schur decomposition of the matrix J gives the vectors qf = 
Vf/\vf \ and qj- = v^/\vf\ directly. The algorithm must be modified to find the 
vectors = Vs/\vs\ and qj- = v^/\vs\, as described in Section |^. | 
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Theorem 4.1 (Planar Case). The equation for the ILDM, Eq. ( \4-di ), admits 
an asymptotic solution in the form of a power series expansion, 

^ = ^(y,5) =^(0)(y)+e^W(y)+e2^(2)(^) + ... as eiO. (4.7) 
The functions ip^^\ and ip^"^^ are defined by the equations 

= ho, (4.8) 
9z^P^'^ = fh'o-9e, (4.9) 

9z 

-\9zM^'^Y-9ze^^'^ -\9ee. (4.10) 

Here, ho = ho{y) is defined by the equation g{y, ho{y),0) = 0, ' denotes dif- 
ferentiation with respect to y, and the functions f and g and their derivatives 
are evaluated at {y,ho{y),0). 



Proof. Assume that z = ilj{y,e), where ip is given by the power series expan- 
sion (|4.7| ). Then 

/(l/,7A(l/,£),£) = / + £(/,^«+/,) 

+ (/.^(^) + yu^^'yf + /.ev^(^) + 1/..) + • • ■ , (4.11) 

where, in the right member, / and its derivatives are evaluated at {y, ip^^^y), 0). 
Similar expansions hold for g and the derivatives of / and g. The leading term 
in the expansion of follows immediately from Eq. ([4.3|) , 

A. = Af + 0(.), Af = /,-^, (4.12) 

9z 

where the derivatives of / and g are similarly evaluated at {y , ip^^\y) , 0) . We 
substitute the various expansions into Eq. ( ^I6| ) and equate the coefficients of 
like powers of e. 

0{e~^). The ILDM equation, Eq. gives 

{gg,){y,z = ^'^^\y),0) = 0, (4.13) 
which is satisfied if tp^'^^ = Hq. This result confirms Eq. {^^1^. 
0{1). From the equation for the ILDM, Eq. ( [4.6|) , we obtain 

fgy + {g,^^'^ +g,)g, = 0. (4.14) 
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Here, we have used the identity g = g{y, hQ{y),0) = 0. The same identity 
imphes that 

gy + gzho = 0, (4.15) 

so Eq. ( [4 .141) reduces to 



0. 



(4.16) 



The assumption of attractive manifolds imphes that gz < 0, so Eq. ( [4.9| ) 
foUows. 



0{e). From the equation for the ILDM, Eq. 



we obtain 



f 



0. 



(4.17) 



Here, we have used Eqs. ( ^4.8| ) and (^4.9|) and the identity ( |4.15|) . The same 
identity also results in a simplification of the expression ( [4.12| ) for A^^-*, 



=fy + f.K. (4.18) 
Furthermore, differentiating Eq. ([4.16| ) with respect to y, we find 

9yzi^^'^ + 9ye + {9zzi^^'^ + 9ze)K + = fK + Uy + fzK)K- (4-19) 

With Eqs. (|4l8| ) and (|l9D, Eq. (|4T7| ) simplifies to 

{9zi''^^ + \9zz{¥'^? + 9ze4^^'^ + \9ee 

- fi^^'^' + {f/9z)K - (/.^^'^ + fe)K) 9z = 0. (4.20) 
Since g^ < 0, Eq. (^A0|) follows. I 



In the following section, we will generalize Theorem 4.1 to the multidi- 



mensional case (Theorem |5.1| ) and compare the asymptotics of the ILDM with 
the asymptotics of the slow manifold (Corollary |0| ). 



5 Asymptotics of the ILDM — General Case 

The definition of the ILDM, Eq. ( |3.4| ), is based on a partition of the Jacobian, 
Eq. (|3.2|) , into a fast and a slow component at each point of phase space 
and a Schur decomposition to generate bases for the corresponding fast and 
slow subspaces. Practical implementations of the Schur decomposition rely 
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typically on the method of deflation [|T^, Chapter 7]; hence, the eigenvalues are 
generated in the order of descending absolute values of their real parts. This 
procedure yields a unitary matrix of the form Q = {Qf Qs)- The columns of 
Q are then reordered, for example by a sequence of Givens rotations, as in 



Ref. 131 



Although this procedure is practical for numerical computations, it is not 
amenable to analysis. We start therefore from the standard Schur decomposi- 
tion before reordering, 

J = QTQ', (5.1) 




^ = h (5-2) 



where 



with Af an n X n upper triangular matrix, an m x m upper triangular 
matrix, and A an nxm full matrix. The diagonal elements of Aj are the 0{e~^) 
eigenvalues of J, and the diagonal elements of A^ are the 0{1) eigenvalues of 
J. The structure of the unitary matrix Q is 

where Qu is an m x n matrix, Q12 an mxm matrix, Q21 an n x n matrix, and 

Q22 an n X m matrix. The columns of ( ) and ( ) form an orthog- 

\ Q21 J \ Q22 J 

onal basis of the fast subspace and its orthogonal complement, respectively. 

Since the fast and slow subspaces are not necessarily mutually orthogonal, 
the orthogonal complement of the fast subspace does not necessarily coincide 
with the slow subspace, and a further operation is needed to identify a basis for 
the slow subspace. This operation consists of solving the Sylvester equation 

AfX - XAs = -A (5.4) 

for the nxm matrix X. With the deflnition 

- ( i ) • (5-5) 

we obtain a block diagonalization of J, 

J = iQY)T,{QY)~\ (5.6) 

where 

T.-('J (5.7) 
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QY 



Qll QiiX + Qi2 
Q21 Q21X + Q22 



-1 



Qll -XQ^2 Q21 -XQ22 
Q12 Q22 



(5i 



Thus, QY reduces the matrix J to its fast and slow components, and the 
condition that the vector field F given in Eq. (|3.1| ) must lie entirely in the 
slow subspace is satisfied if 

(Q'n - XQ[2)f + e-\Q'2, - XQ'22)g = 0. (5.9) 
9D is the same as the ILDM obtained from 



The ILDM obtained from Eq. 
Eq. (|3.4|) and also the same as the ILDMs obtained in Refs. and p4|. 



Theorem 5.1 (General Case). The equation for the ILDM, Eq. i\3.^) , admits 

an asymptotic solution in the form of a power series expansion, 

^ = ^(y,e) =^(o)(y)+e^a)(y)+£2^(2)(^) + ... as e 10. (5.10) 
The -valued functions 'ip^^\ i^^^K and ip^'^'^ are defined by the equations 

= ho, (5.11) 
(D,(?)V'^^) = {Dho)f-ge, (5.12) 
(D,^7)^(^) = iD^P^'^)f~iD^g)-\D'ho)ifJ) + iDho){{DJ)^P^'^ + f,) 

- i(D,2^?)(^«,^«) - {D,gs)^P^'^ - Igse. (5.13) 

Here, Hq = ho{y) is the H"- -valued function defined by Eq. ( l^.^ j, g{y, ho^y), 0) = 
0; DHq = {Dho){y) is a linear operator from R"^ to R", which is represented 
by the nx m matrix of partial derivatives dho^i/dyj, and D'^Hq = D{Dho) = 
{D'^ho){y) is a bilinear map from R™ x R™' to 11^, {D'^ho){u,v) = {{D'^ho)u)v 
for all u,v & R™". The functions f and g and their derivatives are evalu- 
ated at {y, ho{y),0) ; D^f = Dzf{y,ho{y),0) is a linear operator from R" to 
R™, D^g = Dzg{y, hQ{y),0) a linear operator from Ti.^ to R*^, and D^g = 
Dz{Dzg) = Dlg{y, hQ{y), 0) a bilinear map from R" x R" to R". 



Proof. Assume that z = tp{y,e) and that ip is given by the expansion (|5.10|) . 
For the asymptotic analysis of Eq. (|5.9|) , we take 

with Q'^i a unitary n x n matrix and Qn^ an m x n matrix to be determined. 
Thus, Q is unitary to O^e). Higher-order terms can be found in a consistent 
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manner so Q{e) is unitary to any desired order. We take, furthermore, 

= Af(e)=E-'A^f'^+Af^ + --- , (5.15) 

A = A(5) =£-^A(-i)+A(°) + --- , (5.16) 

A, = A,(£) = Af + ■■■ , (5.17) 

X = X(£) = X(o)+£X« + --- . (5.18) 

The generahzation of the expansion (f4.ll ) to the present case is 

f{y, ^{y, e),e) = f + e ((DJ)^« + /,) + e' {{DJ)^^'^ 

+ \{Dlm^'\ ^«) + (D J,)^« + +■■■■ (5.19) 

In the right member, / and its derivatives are evaluated at (y, ■?/'(°^(?/), 0). 
Similar expansions hold for g and the derivatives of / and g. 

To prove the theorem, we substitute the various expansions into Eq. ( |5.9| ) 
and equate the coefficients of like powers in e in the usual manner. 

0{e~^). The ILDM equation, Eq. (^, gives 

Q'-'h = 0. (5.20) 
Since Q^2i is unitary, Eq. ( |5.20|) reduces to 

(7(y,V^W(l/),0) = 0. (5.21) 
This equation is satisfied if = /iq, which confirms Eq. ( ^.11[) . 
(9(1). From the ILDM equation, Eq. ( ^.9|) , we obtain 

- X(o)/ + g(f ((Z^.(7)^(^) + ge) = 0. (5.22) 
Here, we have already used the identity g = 

The matrix X^^^ is determined from the 0{e~^) terms in the Sylvester 
equation, Eq. ( p.4|) , 

^(-1)^(0) ^ (5^23) 

The matrices Aj"^"* and A^^^\ in turn, follow from the 0{e^^) terms in the 
Schur decomposition, Eq. (|5.1| ), 



Q(?A(-i) = Dyg, Q^h^f'^Q^I = D.g. (5.24) 

The second equation is the Schur decomposition of D^g, so Q^2i is determined 
by the ordering of the elements of A^."^''. Both equations can be inverted. 
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A(-i) = Q^^/{Dyg), A[-'^ = Q^hDM^l (5.25) 



Hence, 

= -Q^S^'iD^gy'iDyg). (5.26) 

We can simplify this expression if we use the identity g{y, ho{y), 0) = 0, which 
holds for all y. Upon differentiation, the identity gives a relation between DyQ 
and D^g, 

Dyg + {D,g){Dho)=0. (5.27) 

Note that this is a relation in the space of linear operators from to R". 
With this identity, Eq. (|06D becomes 

xW = Q(°/(mo), (5.28) 

and Eq. (|5.22|) reduces to 

g^f [{D,g)ij^'') +ge- iDho)f] = 0. (5.29) 

Since Q^i is unitary, Eq. ( ^.12|) follows. 

0{e). From the ILDM equation, Eq. (|5.9|), we obtain 

(g(i)' _ x(i) + Qf/iDho)Q^^^QfhDho)) f - Q^hoh,) {{DJ)^^'^ + /.) 

+ Q?/ {{D.gW'^ + UDl9)ii^^'\i^^'^) + iD^9e)i^^'^ + kdee) = 0. (5.30) 

Here, we have already made use of Eqs. (|5.11|) and (|5.12| ) and substituted the 
expression ( p.28| ) for X^^\ 

The matrix X^^^ is determined from the 0{1) terms in the Sylvester equa- 
tion, Eq. (|]4D, 

A^-')X(^) + Af - X(°)Af = -A(o). (5.31) 

The matrices A^''', A^^-*, and A^^-' follow in turn from the 0{1) terms in the 
Schur decomposition, Eq. (|5.1|), 

gS^A(-i) + Af = Dyf, (5.32) 

Q^A^-^^Q^f = DJ, (5.33) 

(A(-^)gff ' + AW) = {D,iDyg))ij^'^ + Dyg,, (5.34) 

gg^ {-A'-'WQf/ + Af^Qf/) = {Dlg)^^^^ + D^g,. (5.35) 

We proceed as follows. First, we solve Eq. ( p. 331 ) for g^ , 

gSl^ = iDJ)Qf^ (A^r'T' = iDJ)iD^g)-'Qfl (5.36) 
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Then, we obtain A^°^ from Eq. ( [5.32| ), 

Af = Dyf - gff A(-i) = Dyf + {DJ){Dho). (5.37) 
(We have used the relation ( ^.271 ) to rewrite the expression ( F25D for A(-i).) 



Next, we solve Eqs. i WM) and (IS^SSp for A^o) and A^^ 



= Q^/ {{D,{Dyg))i^^^^ + Dyg, - (D^g) [{DJ){D,g)-^)'^ , (5.38) 

Af = Q^' {{Dlg)^''' + D.g, + Q^^A^^^^Qff Q?/) Q^^ 

= Q?/ {iDlg)^IJ^'^ + D^g, - {D^g){Dho){DJ){D,g)-^) Qg^(5.39) 

After these steps, we find X*^^) from Eq. ( [5.31[ ), 

X« = (A(-^))-i(-A(°)-Afx(°)+xWAf) 



(D2(?)V'(') + - {D,g){Dh,){DJ){D,g)-^) {Dh 
+ {Dho) {Dyf + iDJ){Dho)) 



(5.40) 



Substituting from Eq. (^36|) and X(^) from Eq. ([CTp into Eq. (|CTD , we 
obtain 

Qf/ [{D.g)-' {mDygm^^'\ f) + {Dyg,)f + {DlgW\ {Dh,)f) 
+ {D,g,){{Dho)f) - {Dh,) {Dyf + {DJ){Dho)) f) 

- {Dh,) ((D,/)^« + /,) 

+ (D.^7)^(2) + \{Dlg){i;^'\ ^(1)) + (D.(?.)^(^) + \gee] = 0. (5.41) 
The bilinear maps Dz{Dyg) and D^^g satisfy the symmetry relations 

{Dz{Dyg)){u,v) = {Dy{D,g)){v,u), n G R^ ^; G R'", (5.42) 
{Dlg){u,v) = {Dlg){v,u), M,t;eR", (5.43) 

so Eq. ( |5.41|) is equivalent with 

Qfl [{D.g)-' {{Dy{D,g)){f,i;^'^) + {Dyg,)f + {Dlg){{Dho)f,i^^'^) 
+ {Dzg,){{Dh,)f) - {Dho) {Dyf + {DJ){Dh,)) f) 

- {Dho) {{DJ)^^'^ + fe) 

+ (Z^.<7)^(2) + \{Dlg){i;^'\i;^'y) + {Dzg,)i;^'^ + \gJ = 0. (5.44) 



We simplify this expression by means of Eq. ( p.l2| ). Upon differentiation, this 
equation gives the identity 

= (D'/io)/ + {Dh^) (Dyf + iDJ){Dho)) . (5.45) 

This is a relation in the space of linear operators from to R". When 
applied to the vector /, it gives the identity 

{Dy{D,g)){f,^|J('^) + {Dyg,)f + {Dlg){{Dh,)f,^(')) + {D,gMDho)f) 
+ iD,g)i{DijW)f) = (D'ho){f, f) + (Dho) {DJ + {DJ){Dho)) f. (5.46) 

With this result, Eq. (|5.44| ) simplifies to 

QS^' [{D.9)i^^'^ + UD^,g){^p('\^P^'^) + {D^g,)ijW + ig^^ 
- (Z^^W)/ + iD,gy\D'ho)if, f) - (Dho) ((/^./)^« + fe)] = 0. (5.47) 

Since Q^^ is unitary, Eq. ( ^.13|) follows. I 



The following corollary summarizes the result of the asymptotic analysis. 

Corollary 5.1 The ILDM is an approximation to the slow manifold Aig of 
the fast-slow system of Eqs. which is asymptotically accurate up 

to and including the 0{e) term as e i 0. The approximation is asymptotically 
accurate up to and including the 0{e'^) term if and only ifD'^ho^y) = for ally. 
The asymptotic expansion of the ILDM is given by Eq. ( ^.IQ ). A comparison 
of the coefficients in the expansion with the coefficients in the expansion of the 
slow manifold M.e, Eq. (\2.12j , shows that 



^'-'^ = ho, (5.48) 
= h^^\ (5.49) 
^(2) ^ h'^^'>-{D,g)-\D'ho){f,f). (5.50) 

The difference ip^'^'^ — h^"^^ involves the bilinear form D'^Hq, which is pro- 
portional to the curvature of the zero set of g a.t e = 0. It is present in 
any fast-slow system, unless the curvature vanishes everywhere. Because of 
it, the ILDM is in general not invariant under the dynamics of the system of 
Eqs. (P)-( 
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6 The Iterative Method of Fraser and Roussel 



The iterative method of Fraser and Roussel was developed originally for planar 
fast-slow systems that are linear in the fast variable, 

y = fi{y,£)z + f2{y,£), (6.1) 

= gi{y,e)z + g2{y,e). (6.2) 

Here, y and z are scalar-valued functions of time. These systems of equations 
are typical for enzyme kinetics |jl2| and other biochemical systems whose dy- 
namics can be reduced to slow manifolds. In this case, the slow manifolds are 
curves in the phase plane. 

On any trajectory z = z{y,e) in the phase plane, we have the identity 
z = Zyij, or, in terms of the functions / and g, 

ezy{fiz + /a) = giz + g2. (6.3) 

This identity is known as the invariance equation. In the present case, it can 
be solved for z in terms of y and Zy, 

z = -S^ + 'J^'y . (6.4) 
gi - ejxZy 

The equation holds, in particular, along trajectories on invariant manifolds. 
Fraser used Eq. ( |6.4| ) to propose the following functional iteration procedure 
to approximate the slow manifold. 

Starting from an initial function one computes a sequence of func- 
tions {(^^ : £ = 1,2,...} using the definitions 

^-JillMtl^ «-1.2..... (6.5) 

9\ - ^fm 

Under appropriate conditions, the sequence {ip^^\y,e) : i = 1,2,...} ap- 
proaches z{y,e) (in a sense to be made precise) as i goes to infinity, so the 
algorithm generates successive approximations to a slow manifold. 



The iterative procedure generalizes to the fast-slow system of Eqs. ( 2.1 



( |2.2| ). The invariance equation is 

£{Dyz){y,e)f{y,z{y,e),e) = g{y, z{y,e),e), (6.6) 
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for any trajectory z = z{y, e) in phase space. Starting from a function one 
computes a sequence of functions {(p^^^ : i = 1,2, . . .} hy solving the equation 

eiDy^^'-'^)iy, e)f{y, ^^'\y, e),e) = g{y, ^^'\y, e),e). (6.7) 

The sequence {(f^^\y, e) : i = 1,2, . . .} approaches z{y, e) (again, in a sense to 
be made precise) as i goes to infinity. 

Notice that Eq. ( |6.7|) amounts to an implicit definition of ip^^\ unless both 
/ and g are linear in the fast variable z, as in the planar case discussed above, 
Eqs. ( |6.1| )-( pl^ ). Hence, the numerical computation of ip^^^ generally requires 
the solution of a nonlinear equation. 



7 Asymptotics of the Iterative Method 

Beacuse the iterative method of Fraser and Roussel is closely related to the in- 
variance equation, its asymptotic properties are most easily analyzed in terms 
of those of Eq. ( |6.6|) . 



Lemma 7.1 The invariance equation, Eq. dj , admits an asymptotic solu- 
tion in the form of a power series expansion, 

z{y,e) = z^''\y)+ez'^^\y) + --- as e i 0, (7.1) 

where 

z^'^ = ho, (7.2) 

and the functions z^''^ : K R", / = 1,2, ... , are found successively from 
Eq. ( \7.1'\ ) below. In particular, z^^^ and z^'^'^ are found from the equations 

{D,g)z^'^ = {Dz^'^)f-ge, (7.3) 
(D,^?)^(2) = (D2«)/ + (Z}^W)((DJ)^« + /,) 

- \{Dlg){z^'\ ^(1)) - {D^g,)z^^^ - \g,,, (7.4) 
where f and g and their derivatives are evaluated at (y , z^^\y) , 0) . 

Proof. We begin by expanding the function /, 

f{y, z{y, e),e) = f^'\y) + ef^'\y) + f^'\y) + ■■■, (7.5) 

where 

f^'\y) = f{y,z^'\y),Q) (7.6) 
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and 

I 



j=l-^' \i\=l 

k=i j=i J- \i\=i~k 
+ l^(dU)(y,z^'\y),0), 1 = 1,2,.... (7.7) 

The derivative {Dif){y, z^'^\y),0) in the first term is a j-hnear map from 
(R")-' to R'". The inner sums are taken over all multiindices i = [ii, . . . , ij) of 
positive integers ii, . . . ,ij with length \i\ = J2i=i ik = I and / — k, respectively. 
The first few coefficients are 

= + (7.8) 

f^\y) = {DJ)z(^^ + \{Dlf){z^'\z^'^) + (DJ,).« + i/.,, (7.9) 
f^'\y) = {DJ)z('^ + {D'j)iz('\z('y) + l{Dlf)iz^'\z^'\z^'^) 
+ {D.Mz^'^ + \{D%){z^'\z^'^) 

+ i(D,/,,)^« + (7.10) 

where / and its derivatives are evaluated at {y, z'^'^\y),Qi), and the argument 
of each z^''\i = 1, 2, . . . , is ?/. A similar expansion holds for g{y, z{y, s),e), 

g{y, z{y, e),e) = g'^^\y) + eg'^'\y) + e^g^^\y) + ■■■, (7.11) 

where 

g^''\y)=g{y,z^'\y),0), (7.12) 

and 

g^Hy) = j:\{Dlg){y,z^'\y),0)Y.{^^^'\y).---.z^^^\y)) 
j=i^- \i\=i 

+ E^E:!r(^^(^'^?))(2/'^^°H2/)'0) E {z^''\yl---,z^'Hy)) 

k=l j=l J- \i\=l~k 

+ ^{dig){y,z'^'\y),0), 1 = 1,2,.... (7.13) 

Termwise differentiation of the asymptotic expansion ( [7.1|) gives 

{Dyz){y, e) = Dz^^^ + eDz^^~^ + ■■■ . (7.14) 
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Equating the coefficients of like powers of e in the left and right members of the 
invariance equation, Eq. ( |6.6|) , we obtain a sequence of functional identities, 



9^'^ = 0, (7.15) 
= Y.iDz^'^)f^'~'"\ £=1,2,.... (7.16) 

We satisfy the 0{1) equation, Eq. ( [r.l5| ), by taking z^^^ = ho; see Eq. (fO]). 
Then f^^\ y) = f{y,ho{y),0), and f^^\y) and g'^'^\y) are given by Eqs. (|77|) 
and ( |7.13| ), respectively, with z^^\y) replaced by ho{y). 

Next, we turn to the 0{e^) equation, Eq. (|7.16| ). We observe that z^^^ 
occurs in g^^^ only in the ffist term with j = 1; the remaining terms involve 
z^^'^ through but not z^^\ The right member of Eq. ( [7.16|) similarly 

involves z^^^ through z^^^^^ but not z^^\ Therefore, the identities (|7.16|) can 
be solved successively for and so on. Thus we find 

- E ^ E UdM9)) E (^^"^ • • • ' ^^''^) - Ws9) (7.17) 
fc=i j=i J- \i\=e-k 

for i = 1,2,.... Here, the functions / and g and their derivatives are evaluated 
at {y, z = ho{y), 0), and it is understood that a sum is empty when the lower 
bound exceeds the upper bound. The equations for £ = 1 and i = 2 are given 
in the statement of the theorem. | 



The following theorem shows that i successive applications of the iterative 
algorithm of Eraser and Roussel, starting from Lp^^^ = ho, generate an approx- 
imation v?*^^^ to the slow manifold Aie that is accurate up to and including the 
0{e^) term as e | 0. 

Theorem 7.1 Let (p^^^ and z^^^ be defined recursively for i = 1,2,... by 
Eq. ( jgT^ j and \7.11[ ), respectively. If (f^^'' = 2^°^ = ho, then 

^ y.W(l/,5) = j:e^z^Hy) + Oie'^'), £ = 1,2, ... . (7.18) 

i=0 

Proof. The proof is by induction. Taking £ = 1, we have 

eiD^^'y)iy)fiy, ^«(y, e),e) = g{y, y.«(y, e),e). 
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Since y?*-''-' = z^^\ this equation is the same as 

eiDz^'^)iy)fiy, e),e) = g{y, ^^'\y, e),e). 

We expand the terms in this equation in powers of e and equate the coefficients 
of hke powers of e. To leading order, we find the equation 

g{y,V^'\y,0),0) = 0, 

which is precisely the equation for z^^\ so 

^«(l/,0) = zW(y). 

To the next order, we find the equation 

iD,g)iy, z^'\ 0)4'\y, 0) + g^iy, z^'\0) = {Dz^'^){y)f{y, z^'\y), 0), 

which is precisely Eq. { \l.3\) , so 

^«(l/,0) = z«(y). 

Thus, 

y.«(l/,e) = zW(2/) + £z«(2/) + 0(£2), 
and the theorem is true for i = 1. 

Suppose the theorem is true for £ — I, y^^^^-^^ = J2i=le^z^^^ + 0{e^). The 
function (^W is defined by Eq. ( |6.7|) , 

eiDy^^'''^)iy, e)fiy, ^^'\y, e),e)= g{y, ^^'\y, e),e). 

We expand each term in powers of e and equate the coefficients of like powers, 

J2{DyZ^'^)f^'-'-'^ =9^^\ j = l,2,...,£, 

where /^'^ and g^'^ are the functions defined after Eq. (fTS]). For j = 1, . . . , i, 
these are exactly the same functional identities as we derived above. Hence, 
order by order, the solution is 

9:^w(y,o) = ^«(y), 2 = 0,1,. 

This proves that the theorem is true for i. | 

Remark. In general, the coefficient of in (p^^^ will not be equal to z^'^\ This 
may be seen by direct examination of the equation, 

+ fs) = \{Dlg){z^'\ z^^^) + {D^gM'^ + \g,,, 
24 



which differs from Eq. ( [r.4| ). The latter has two additional terms, namely, 
{Dzg)z^'^^ and {Dz^^^)f. Hence, the first iterate generally involves an error of 
0{e'^). Similarly, the ith iterate does not give the same equation for the term at 
order 0{e^~^^) as compared with that obtained from invariant manifold theory; 
hence, the error in the approximation at this stage is generally 0{e^~^^). I 



A comparison of the results given in Theorem |7]1| with Theorem |2.1| leads 
to the following conclusions. 

Corollary 7.1 The iterative method of Fraser and Roussel gives successively 
higher-order asymptotic approximations to the slow manifold Ai^. Starting 
from Lf^^^ = Hq, i applications of the iterative procedure give an approximation 
Lf^^^ that satisfies 

= ^£*/i» + C(£^+i), £ = 1,2,..., (7.19) 

1=0 

where the functions /i*^*^ are the coefficients in the asymptotic expansion of M.^ 
given in Theorem 2J_, Eq. ^2.12J . 



Remark. In Ref. , Roussel and Fraser noted the connection between 



the iterative approach developed by Fraser and Roussel and techniques from 
dynamical systems theory, such as deriving the existence of certain invariant 
manifolds via application of the contraction mapping principle to appropriate 
integral equations. I 



8 The Michaelis—Menten— Henri Model 



In this section, we illustrate the analytical results of the preceding sections on 
the Michaelis-Menten-Henri (MMH) model. The MMH model is a prototype 
reaction mechanism for enzyme kinetics in biochemistry [^. It is a planar 



fast-slow system, which is given in nondimensional form by the equations 

y = -y+{y + a-b)z, (8.1) 
ez = y-{y + a)z. (8.2) 

The variables are y, the concentration of the substrate, and z, the concentra- 
tion of an intermediate substrate-enzyme complex. The parameters a and b 
satisfy the inequalities a > 6 > 0. 
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Remark. The MMH model, Eqs. ( ^.I| )-( pl2[ ), is derived from a more com- 
plicated system involving four species (enzyme, substrate, enzyme-substrate 
complex, and product) and two reactions (one reversible, one irreversible). The 
full system can be reduced to the planar system, because it has two conserved 



quantities; see Ref. MM. | 



The system of Eqs. (p.l|)-( p^ has a family of slow manifolds Ad^, whose 
asymptotics are given by Eq. ( p.l2|) , 

h,iy) = hoiy)+eh^'\y)+e'h^'\y) + --- , y > 0, (8.3) 

where 

hoiy) - (8.4) 

'''^y^ - (^' ^^•^) 

h(^\y) = -byi2ab-3by-ay-a^)_ 
{y + ay 

We now show that the ILDM method yields an approximation of the slow 
manifold that is accurate up to and including the 0{e) term, but differs at 
0{e^) because of the curvature of Hq. 

The Jacobian of the vector field associated with Eqs. ( ^.l|) -( p^ is 

T_ / -(1-2;) y + a-b \ . . 

and its eigenvalues are 



The (nonnormalized) slow eigenvector is 



Xs + £ ^{y + a) 



.9) 



and there is a corresponding fast eigenvector. The vector Vg spans the slow 
subspace. As noted before, the fast and slow subspace are not orthogonal, so 
we work with vj- and define the ILDM as the set of points where the vector 
field is orthogonal to vj-, 

(1 -z)[-y+iy + a- b)z] - (A, + e-\y + a))[y - (y + a)z] = 0. (8.10) 
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Asymptotically, the ILDM is given by 

z = ij{y, e) = V^W(y) + 62P^'\y) + e'ij^^\y) + ■■■ , (8.11) 



where 



i^^'\y) = (8.12) 

y + a 



(y + a 



i.i^\y) = ^byi2ab-by-ay-a^)_ 

A comparison of the coefficients in the expansions ( ^.3[ ) and ( |8.11| ) shows agree- 
ment of the 0{1) and 0{e) terms. On the other hand, the 0(£^) terms differ; 
their difference is proportional to the curvature of ho, 

^(^) - /.(^) = = (8.15) 

[y + ay g; 

The iterative method of Fraser and Roussel starts from the invariance equation, 

y + £Zyy 



y + a + ezy{y + a — b) 



(8.16) 



or, equivalent ly. 



2 



y , byzy by{y + a-b)z 



3 



y + a (y + ay {y + a) 

o byiy + a-bYzl 

+ ^7 xF^^T T^- 8.17 

{:y + aY + e{:y + aY{y + a-b)zy 

Successive applications of the iterative algorithm lead to the approximations 

(^(°) = (8.18) 
y + a 

^^'^ = ^ + e-^, - ,^-^Mi±^ + o{^), (8.19) 
y + a {y + a)^ {y + ay 

^(2) = -by ^ ^,abyi2ab-3by- ay-a^) ^ ^^ ^0) 

y + a {y + ay [V + ay 

A comparison with Theorem |2.1| (and Eqs. ( ^.4|) -( ^^ ) shows that is 
asymptotically correct up to and including terms of 0{e^) for £ = 1,2 (and 
beyond), as predicted by the analysis. 

Remark. The analysis of the full MMH model has been significantly extended 
in Ref. | 
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9 The Davis-Skodje Model 



The planar fast-slow system 

y = -v. (9.1) 

ez = -^ + 7^-77^, (9.2) 
1 + y [l + yr 

was introduced by Davis and Skodje as a model on which to compare 
various reduction methods. (The inverse e^^, which is large, equals the large 
parameter 7 of Ref. 0.) 

For any e, the curve 

z = hM = -r^, y>0, (9.3) 
1 + y 

is invariant under the dynamics of Eqs. ( p.l| )-( p^ ). Therefore, the function 
hs represents the slow manifold exactly on y > for all small e > 0. (The 
nonlinearity in Eqs. ( |9.1| )-( PT^ ) was, in fact, chosen so the slow manifold is 
given by the simple expression of Eq. ( p.3|) .) 

The Jacobian of the vector field is 

'^^{e-\{l + y)+e{y-mi+y)-' -^"0 ' ^^'^^ 

and the eigenvalues are = — 1 and A/ = —1/e for all {y, z). The correspond- 
ing (nonnormalized) eigenvectors are 

{I - e)-\l + y + e{y - + y)-^) ^ 
The ILDM is given by the expression 



^/ = 1 • (9.5) 



z = -^{y, e) = —— + ; 9.6 

l + y {l-e){l + yy 



cf. Eq. (3.8)]. Its asymptotic expansion is 



^ +^'77^ + ---- (9.7) 



l + y (l+yY 

The error in the expansion is 0{e'^) and proportional to the curvature, h'^ 
-2/{l+yf. 
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The invariance equation is 



" = TT^-^(rT^ + ^^"- (9-8) 

from which one readily verifies that the iterative method of Fraser and Roussel 
yields the approximation ^^^\y) = y/{l + y) for £ = 1, 2, ... . 

Remark. If one restricts the variable y to a finite interval, then the sys- 
tem of Eqs. (|9.1|) - (|9.2|) has a family of slow manifolds, all exponentially close 
(0(e"'^/^^) for some c > 0) to the exact slow manifold, Eq. (|9.3| ). | 



10 Discussion 



The fast-slow system of Eqs. (|2.1|) - (|2.2|) captures the essential elements of any 
reaction mechanism whose long-time dynamics evolve on a slow manifold in 
the composition space. As stated in Section |^, however, this system is a math- 
ematical idealization, and we need to consider how the results of the analysis 
carry over to more general reaction mechanisms. In this subsection, we con- 
sider issues related to the separation of time scales (Section p.0.1| ), the inclusion 
of conserved quantities (Section |10.2|) , and the development and analysis of re- 
duction mechanisms for reaction-diffusion equations (Section p.0.tj| ). 



10.1 Separation of Time Scales 

In this section, we discuss the partition of variables into a fast and a slow 
group and the assumption that the groups evolve on time scales that are and 
remain well separated at all times. This assumption underlies the definition 



of the small parameter e in the model of Eqs. ( |2.1| )-( |2l^ ). We also discuss the 



possibility of partitioning the variables in more than two groups. 

While many of the systems to which reduction methods have been applied 
satisfy this assumption, there are a significant number of reaction mechanisms 
where the fast and slow time scales are separated, but not well separated. 
This is the case, for example, when e is no longer an asymptotically small 
parameter but a fixed (relatively small) number. In such cases, the spectral 
gap between the fast and slow eigenvalues of the Jacobian of the vector field is 
small, much smaller than the chasm between the fast 0{6~^) eigenvalues and 
the slow 0{1) eigenvalues for asymptotically small e, even at points near a low- 
dimensional manifold. Nevertheless, as long as there exists a nonzero spectral 



29 



gap, the Jacobian can still be reduced to a fast and a slow component, and 
the ILDM method can be (and has been) implemented numerically. Part of 
our ongoing research is aimed at using spectral projection operators to analyze 
these applications of the ILDM method. 

In addition, there are reaction mechanisms where the number of fast and 
slow species changes over time, as may happen, for example, when the tem- 
perature in the chemical reactor changes and the least-slow of the slow species 
transits to the group of fast species. To account for this type of occurrence, 
practical implementations of the ILDM method may use one ILDM until the 
crossover occurs and another after the crossover. The two procedures can be 
linked by doing a numerical integration of the full system of equations during 
the crossover. The analysis presented in this article applies to each ILDM 
separately. However, it should be noted that the relationship between these 
two ILDMs (and, more generally, slow manifolds) depends on the bifurcation 
that occurs at crossover. As an alternative, one may consider repartitioning 
the species to avoid crossover altogether. 

In some systems, a component of y may evolve on an even slower time 
scale than that given by r = et. For example, in Eq. ( |2.1| ) one may have 
y[ = e^^^ fi(y, z,e) with 7 > for some index i. Such cases are accounted for 
by the present analysis; in fact, it suffices to absorb the factor e"^ in /j. On the 
other hand, the fact that there are slower time scales in the model may point 
to the existence of still lower- dimensional slow manifolds. By systematically 
eliminating fast variables, starting with the fastest and proceeding up the hier- 
archy, one can reduce the dimensionality of the slow manifold in a systematic 
way until no further reduction is possible. Such an approach has been used, 
for example, in Refs. 0, |7^. The idea of a hierarchy of time scales was first 
explored in singular perturbation theory by Tikhonov [K9 . 



10.2 Conserved Quantities and the ILDM Method 

In this section, we briefly consider the ILDM method for fast-slow systems 
whose dynamics are described by Eqs. (|2.1| ) -(|2.2| ) , where the unknowns satisfy 
one or more conservation laws. For models of chemical reactions, a conserved 
quantity is, typically, a linear combination of several unknowns that is constant 
in time. Conserved quantities give rise to zero eigenvalues of the Jacobian of 
the vector field, and the ILDM method groups these zero eigenvalues with the 
slow ones. 

Degeneracies in the Jacobian affect the analysis of the ILDM method. 



Consider, for example, a system given by Eqs. (|2.1|) - (|2.2|) with m = 2 and 
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n = 2, which has one conserved quantity. If the conserved quantity is a hnear 
combination of the two slow variables, then the first and second row of the 
Jacobian are linearly dependent. In this case, there is effectively only one 
slow variable, and the analysis presented in this article applies to the one- 
dimensional slow manifold Aio. If, by contrast, the conserved quantity is a 
linear combination of the two fast variables, then not only is J degenerate, 
but also the columns of Dzg{y,ho{y),0) are linearly dependent; hence, A4q 
is not asymptotically stable in the four- dimensional composition space. A 
reduction of the number of fast variables by one lifts this degeneracy, because 
A^o is asymptotically stable in the three-dimensional reduced composition 
space. Other possibilities are that a conserved quantity depends on a mix of 
fast and slow variables or that there are multiple conserved quantities. In each 
case, the analysis of the ILDM method presented here applies once the system 
has been reduced to a system for which the manifold A^o is asymptotically 
stable. 



10.3 React ion- Diffusion Equations 



Reduction methods have also been developed for systems of reaction- diffusion 
equations [|^, |7^, ^ [19|, |6^, |3^. The elimination of the fast species 
affects not only the reaction kinetics of the slow species (as is the case for 
the kinetics reduction methods considered in this article) but also their dif- 
fusivities. The "effective" diffusivities differ from the ordinary diffusivities by 
concentration-dependent terms that are higher order in e. 

We illustrate this phenomenon on the Michaelis-Menten-Henri mecha- 
nism with diffusion of the slow species. 



y = -y+{y + a-h)z + DAy, 
ez = y-{y + a)z. 



(10.1) 
(10.2) 



The variables y and z depend not only on time but also on space; A is the 
Laplace operator. Ideas from inertial-manifold theory have been applied to 



this react ion- diffusion system |7^. The slow manifold is infinite dimensional, 
its asymptotics are given by 



y 



+ e 



aby 



:DAy 



+ 0{e' 



(10.3) 



y + a ' ' [(?/ + a)4 {y + a)^' 

and the reduced react ion- diffusion equation up to and including terms of 0{e) 
is 

^^y \ : f ^ ^(y + 0, — b) 



y = -y+{y+a-b) 



y I ^ 

y + a ""{y + a) 



Dll-e- 



(y + a) 



Ay. (10.4) 
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The regular diffusivity D is seen to be corrected by an e-dependent term that 
involves the concentration of the slow species. 

Singular perturbation theory provides an alternative method to find the 
infinite-dimensional slow manifold and, hence, the reduced reaction-diffusion 
equation for this and similar systems. In particular, one finds an asymptotic 
expansion for the slow manifold of the form 

z = ho{y) + eh^'\y, Ay) + e^h^^\y, Ay, A^y) + ■ ■ ■ . (10.5) 

For the MMH model, this procedure leads to the same expansion of the slow 
manifold, Eq. ( |10.3| ), and dynamical systems theory states that is the 
infinite-dimensional weak stable manifold of the spatially homogeneous state 
(0,0). 

We remark that one can include diffusion of the fast variable z, but only 
if the diffusion coefficient is 0{e), 

y = fiy,z,e) + D,Ay, (10.6) 
ez = g{y,z,e) + eD2Az. (10.7) 

In this case, z^^^ is independent of D2, and one can follow the same asymptotic 
procedure as before; the diffusion coefficient D2 enters into the equation of 
order 0{e). However, if the diffusion term in the fast equation is 0{1), the 
asymptotic procedure for finding z^^^ breaks down. 
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